Bees in the six: Determinants of bumblebee habitat quality in urban landscapes

Abstract With growing urbanization, it is becoming increasingly important to design cities in a manner that sustains and enhances biodiversity and ecosystem services. Native bees are critical pollinators that have experienced substantive declines over the past several decades. These declines have captured the attention of the public, particularly urbanites, prompting a large interest in protecting pollinators and their habitats in cities across North America and Europe. Unfortunately, we currently lack research about specific features of urban environments that can enhance the fitness of pollinators. We carried out an intensive study of Bombus impatiens, the Common Eastern Bumblebee, in the city of Toronto (Canada's largest city), to better understand landscape parameters that provide high‐quality habitat for this species and likely other generalist bees. We divided the city into 270 grid cells and sampled a large number of worker bees, which were then genotyped at twelve hypervariable microsatellite loci. The genetic data allowed us to quantify the effective number of colonies and foraging distance for bumblebees in our study area. We then asked how the city's landscape and human population demography and income are associated with the availability of high‐quality habitat for B. impatiens. Several aspects of Toronto's landscape influenced colony density and foraging range. Urbanization had a clear effect on both colony density and foraging distance of workers. On the other hand, functional (i.e., not cosmetic) green space was often associated with higher quality habitats for bumblebees. Our study suggests several planning strategies to enhance habitat quality for bumblebees and other pollinators in cities.


| INTRODUC TI ON
Approximately 55% of the world's human population live in urban areas, and this is expected to grow to nearly 70% by 2050 (UN, 2018). Urban environments can sustain a large number of native, and sometimes at-risk, species thereby providing considerable value for conserving biodiversity and ecosystem services (Nilon et al., 2017).
Moreover, cities contain a substantive proportion of the world's human population that-if exposed to nature and wildlife-can act as catalysts and ambassadors for species conservation. However, it is not always obvious how to design cities in a manner that best supports biodiversity and ecosystem services (Lepczyk et al., 2017).
There has been a growing concern over the decline of bees across the globe (Cameron et al., 2011;Grixti et al., 2009;Potts et al., 2010;Powney et al., 2019;Senapathi et al., 2015). These declines have captured the attention of the public, leading to widespread interest in conserving and protecting bee populations. The free ecosystem services provided by wild bees in cities are important for continued pollination of community vegetable gardens, native plant gardens, indigenous medicine gardens, naturalized green spaces, urban agricultural fields, and fruit trees. Cities can host diverse pollinator assemblages (Hall et al., 2017), and several major North American and European cities have implemented mandates or legislation to protect pollinators. Indeed, cities can now seek pollinator-friendly certifications through programs such as "Bee City" (e.g., beecitycanada.org and beecityusa.org) by pledging to implement a suite of actions aimed at protecting pollinators. While some of the general recommendations for protecting pollinators in cities make intuitive sense, they are not often backed by research. Planting "pollinator gardens" is often the top recommendation but we note that its efficacy in improving bee health is rarely empirically demonstrated (but see Baldock et al., 2019;Leve et al., 2019). Installing "bee hotels" (i.e., fixtures with hollow tubes for bees to nest in) was a common recommendation for conserving pollinators in cities but it has since been shown to not benefit native bees (MacIvor & Packer, 2015). We need research to understand which features of urban environments enhance pollinator health to help inform bee conservation initiatives in cities.
Several studies in mixed environments have identified the importance of landscape as a driver for bee diversity and habitat quality. For example, sampling twelve sites in urban, seminatural, and agriculture areas in France, Geslin et al. (2016), discovered that "urbanization"-quantified as proportion of impervious surfaces within a site-was negatively associated with bee abundance and species richness. Similar patterns were observed within farms and gardens in the greater Chicago metropolitan area (Bennett & Lovell, 2019). Using genetics to infer the location and density of bumblebee colonies, Goulson et al. (2010) discovered that the size and number of gardens within agricultural landscapes had positive effects on the survival of two bumblebee species in the UK. Another UK study carried out in agricultural landscapes found that forage patches sown with flower mixtures had a positive influence on bumblebee species richness and worker densities (Carvell et al., 2011). Jha and Kremen (2013a) found that floral diversity across mixed landscapes (i.e., crops, grasslands, orchards, riparian forests) was associated with the foraging distance for the bumblebee Bombus vosnesenskii.
In the same study, the authors found that the proportion of paved surfaces had a negative influence on Bombus vosnesenskii's nest density. The above studies in mixed or agricultural landscapes suggest bee habitat quality is enhanced by the availability of floral resources and is negatively impacted by impervious paved surfaces.
However, it still remains to be demonstrated if the same drivers of habitat quality identified in studies of mixed or agriculture landscapes still operate at the within-city scale that is needed to inform urban planning.
In addition to physical features, the characteristics of the human populations inhabiting urban landscapes can potentially influence the quality of bee habitats via a phenomenon known as the "luxury effect" (Leong et al., 2018); wealthier neighborhoods may actively or passively provide higher quality habitats for wildlife than poorer neighborhoods. A recent study discovered that higher household income was positively associated with pollinator abundance in the UK (Baldock et al., 2019), although household income was not associated with bee abundance, richness, or community composition in Chicago (Lowenstein et al., 2014). Understanding the importance of human population demography more broadly in pollinator conservation will be critical for ensuring that ecosystem service provisions in cities are equitable.
Here, we leveraged Toronto's (Ontario, Canada) rich open data resources to study how the physical features of the city and the demography of its human population influence habitat quality for bumblebees. Toronto is the capital city of the province Ontario and, with an estimated population size of over 6 million people in the city and its surrounding suburbs (Statistics Canada, 2020), is the largest city in Canada and among the top ten largest cities in North America. We studied the Common Eastern Bumblebee Bombus impatiens (Williams et al., 2014) as it has many traits that may allow us to generalize our findings to other bumblebees and native bees (also see the Discussion section): (1) B. impatiens (Figure 1) is common F I G U R E 1 A common eastern bumblebee, Bombus impatiens, worker foraging on purple coneflower (Echinacea sp.). Photograph by Amro Zayed throughout the city, thereby allowing us to explore the relation between landscape features and bee habitat quality across a large swathe of Toronto; (2) B. impatiens is polylectic (i.e., pollen generalists) like the majority of bees found in cities (MacIvor et al., 2014;Matteson et al., 2008); (3) B. impatiens visits many plant species that are known to be attractive to a wide range of native bee species (Colla & Dumesh, 2010;Gervais et al., 2020;Vaudo et al., 2014); and (4) B. impatiens is social with adults active from early spring to late fall (Colla & Dumesh, 2010), thereby increasing the temporal scope of our analysis. While B. impatiens is certainly one of the largest bees found within Toronto, the reasons described above still suggest that it can act as a proxy for a large number of generalist bees inhabiting the city.
To study how Toronto's landscape influences B. impatiens, we first divided the city into 270 2 × 2 km grid cells and then randomly sampled B. impatiens workers from across the city. The sampled workers were genotyped allowing us to estimate the number of colonies, and how far workers travel to forage, across the different grid cells. We then carried out analyses to examine which physical and demographic features were associated with habitat quality for B. impatiens. In lieu of direct data on reproductive fitness, we indirectly gauged high-quality habitat as features that reduce foraging distances and increase colony density within grid cells. Foraging is metabolically costly which can directly influence the lifespan of workers and fitness of colonies (Kelemne et al., 2019;Rueppell et al., 2007;Tomlinson et al., 2017). Areas with ample floral resources nearby should thus lessen the metabolic costs of foraging allowing bumblebee colonies to allocate more resources toward reproduction, thereby increasing colony fitness. Similarly, we expected areas with more floral resources in urban environments to sustain more bumblebee colonies (i.e., greater colony density per unit area), as previously found in agricultural and mixed-use landscapes (Osborne et al., 2008).

| Population sampling
The city of Toronto (Ontario, Canada) was divided into 270 grid cells, each measuring 2 × 2 km (4 km 2 ) ( Figure 2). Grid cells were selected at random and visited for sampling between July and October 2016this period is within B. impatiens' worker activity period in Southern Ontario (Colla & Dumesh, 2010). In total, we collected 760 specimens from 86 grid cells, across 27 collection days (Data S1). At the start of a collecting day, a random number generator was used to pick a focal grid cell for sampling. We typically traveled within a grid cell and netted flying bumblebees as they foraged on flowers. While we tried to sample bees from many different locations within a grid The exact sampling location where each individual was caught was recorded for downstream analysis.

| DNA extraction
We extracted DNA from bee tissues using the Mag-Bind ® Blood & Tissue DNA HDQ 96 Kit (Omega Bio-tek Inc.) optimized for the KingFisher™ Flex Purification System (Thermo Fisher Scientific Inc.). Briefly, one half of a bee's thorax was finely ground in a 1.5 ml microcentrifuge tube in liquid nitrogen. Tissue lysis buffer (350 µl; Omega Bio-tek Inc.) and Proteinase K (20 µl; Omega Bio-tek Inc.) were added, and the tube was vortexed for 10 s, followed by an overnight incubation at 50°C. Samples were then centrifuged at 7000 g/rcf for 10 min. The supernatant (approx. 300 µl) was transferred to a new tube and centrifuged similarly for another 10 min.
We then followed the kit's published protocol for extracting DNA using the KingFisher™ Flex Purification System. The eluted purified DNA (100 µl) was used for microsatellite genotyping.

| Microsatellite genotyping and scoring
We used fluorescently labeled primers (Table 1) to genotype 12 hypervariable microsatellite loci previously identified for bumblebees (Estoup et al., 1995(Estoup et al., , 1996Funk et al., 2006). Each locus was tested for usability in a two-step process on a small number of B. impatiens workers. First, a temperature gradient was run, with the polymerase chain reaction (PCR) conditions described below, to determine optimal annealing temperatures and product sizes. The gradients were set to 3-5°C above and below the temperatures listed in Estoup et al., 1995, Estoup et al., 1996, and Funk et al., 2006. PCR products were run on 2% agarose gels and genotyped (see below).
Agarose gel results and chromatograms were used to group primers for multiplexing (i.e., multiple loci amplified in a single PCR reaction) or poolplexing (i.e., loci amplified separately and combined for genotyping) based on annealing temperatures, product sizes, and product intensities. Next, the forward primers in each group were labeled with a different color fluorescent tag and each group of multiplexed or poolplexed loci was tested for usability via PCR and genotyping.
After optimization and testing, the 12 primers were amplified as 3 multiplex sets (each with 2 primers) and 2 poolplex sets (each with 3 primers) as shown in Table 1.
All PCRs were performed using a Mastercycler (Eppendorf, Canada) in either strip tubes or 96-well plates with one negative control, where sterile water was used instead of DNA. Each 25 µl PCR consisted of 0.5 µl of each 10 µM primer (Integrated DNA Technologies Inc.), 12.5 µl Taq 2X Master Mix (New England Biolabs Lt.), 1.5 µl DNA, and nuclease-free water (Thermo Fisher Scientific Inc.). PCRs were carried out with an initial denaturation step of 94°C for 3 min; followed by 35 cycles of 94°C for 30 s, T m (Table 1) for 30 s and 72°C for 30 s, with a final incubation at 72°C for 10 min. Following multiplexing and poolplexing, 10 µl samples were sent to The Center for Applied Genomics for genotyping using an ABI3730xL DNA Analyzer (Applied Biosystems) with GeneScan 500 LIZ Size Standard (Applied Biosystems). Geneious v11.1.2 (Biomatters Limited) was used to view chromatograms and to automatically bin alleles and assign genotypes. All assigned bins and genotypes were manually checked for accuracy and adjusted where needed. Micro-checker v2.2.3 (van Oosterhout et al., 2004) was used to identify potential sources of genotyping error. All loci passed Micro-checker's quality control steps with no evidence of scoring error due to stuttering, large allele drop-out, or null alleles.

| Population genetic analyses
We used Colony v2.0.6.4 (Jones & Wang, 2010) to partition worker bees into full sister groups using the following parameters: (1) haplodiploidy; (2) monogyny and single queen mating; and (3) 0% allele dropout rate and 1% rate for other genotyping errors, including mutations (based on our microchecker results). We ran this analysis 3 times using long runs to ensure convergence of the sibship reconstructions. Following established methods (Carvell et al., 2017;Dreier et al., 2014), we retained sibship clusters with probabilities greater than or equal to 80%; this resulted in 303 worker bees grouped into 118 clusters, each containing 2 or more sisters, and 457 singletons (i.e., worker bees that did not have any sisters in our dataset).
We used this dataset to estimate the effective number of colonies and average foraging distance for each grid cell following established methods (Carvell et al., 2017;Dreier et al., 2014;Redhead et al., 2016), as described below. Here, we decided to exclude all 457 singletons from further analysis because we cannot be certain if they resided in a grid cell, or were simply foraging in a grid cell. Moreover, of the 118 sibship groups detected, 31 (totaling 79 worker bees) contained only sister bees with the same geographic coordinates. These groups were also removed from further analysis because they could not be used to calculate foraging distance using our method of triangulation (see below). The remaining 87 sibship groups (totaling 224 worker bees), containing sister bees having 2 or more different (i.e., not identical) coordinates, were mapped in ArcGIS v10.6 (Esri, USA).
Following established methods (Carvell et al., 2017;Dreier et al., 2014;Redhead et al., 2016), colony locations were calculated as the mean center of the sampling coordinates for all sister bees belonging to a sibship group, and the total number of colonies were tallied for each grid cell. However, the total number of colonies we detected for each grid cell represents a minimum. Therefore, we followed Charman et al. (2010) to calculate the effective number of colonies (colNe, i.e., minimum number of sampled + unsam-  2010), colNe in each grid cell was calculated as 1.5 × the total number of colonies in that grid cell, assuming monogyny and monoandry. The first assumption is clearly supported as polygyny is extremely rare in bumblebees (Cameron & Jost, 1998), and the available genetic data on B. impatiens have not detected nests with multiple queens (Payne et al., 2003). While multiple mating is known to occur in B. impatiens, it is relatively rare and the average effective number of mates per queen (1.06) is close to the value expected with monogamy (1) (Payne et al., 2003). Rare cases of multiple mating would lead to undetected halfsib relationships, but we do not expect this to lead to confounds in the downstream analysis. To calculate average foraging distance (aveMeanFD), colony-specific worker foraging distances were calculated as the mean of all the Geodesic (straight-line) distances to a colony for all sister bees belonging to that colony (even if the sisters were captured in different grid cells). Then, grid cell values were calculated as the average of all the foraging distances for colonies located in a grid cell.
We used PopGenReports (Adamack & Gruber, 2014) to estimate basic population genetic parameters (e.g., number of allele, observed and expected heterozygosity) and to test for deviations from Hardy-Weinberg equilibrium at each locus. For this analysis, our dataset was pruned by randomly retaining one sister from each colony that contained 2 or more sister bees. We used a Bonferroni correction when assessing significant departures from Hardy-Weinberg equilibrium.

| GIS and City of Toronto landscape layers
Twenty-eight physical and demographic variables were extracted from 10 maps obtained from the following sources: (1)

| Statistical analyses
We were able to estimate the effective number of colonies and average foraging distances within 49 2 × 2 km grid cells in Toronto (Data S1). We first examined the distribution of the population genetic and landscape parameters and found that transforming the average foraging distance using natural log improved normality (Shapiro-Wilk normality test, p = .07). To understand the relationships between Toronto's landscape and demographic variables (Table 2) and effective number of B. impatiens colonies and average log foraging distance, we performed Spearman correlations using the library Hmisc (Harell, 2021) and corrplot (Wei & Simko, 2017) in R. We corrected p-values for multiple testing using the Benjamini-Hochberg false discovery rate (FDR) adjustment (Benjamini & Hochberg, 1995). We report adjusted p-values in the results section.
We also performed multiple linear regressions using statsmodels (Seabold & Perktold, 2010) to understand how landscape and demographic features of Toronto influence bumblebee effective nest density and average log foraging distances. We first selected the models with the highest adjusted R 2 and then selected the most parsimonious model using the Akaike information criterion, AIC (Burnham & Anderson, 2002), where the lowest value is the most parsimonious model. For the most parsimonious models, we then used the function varpart from the vegan package (Oksanen et al., 2019) to partition the variation between the different explanatory variables.
We then performed a redundancy analysis (RDA) to explore the association between log-transformed foraging distance and number of unique bumblebee colonies with landscape and human demographic variables, using the "vegan" package in R (Oksanen et al., 2019). RDA is a direct gradient ordination approach that allows us to summarize the linear relationship between our explanatory variables and a multivariate set of response variables (Rao, 1964). To account for multicollinearity between explanatory variables, we calculated correlations between variables and used the variance inflation factor (VIF>5) (Dormann et al., 2013) to remove variables that were highly correlated to others. A number of multicollinear variables (tree canopy cover, % city parks, % forests in watersheds, % other paved surfaces, total human population, human population over 60 years old and human population between 40 and 59 years old) were retained in our final model as they were biologically relevant to our study and cumulatively explained significant amounts of variation. We decided to retain these multicollinear variables because removing them would make the biological interpretation and actionable implications of our analysis more difficult. We used a permutation test to determine global significance of the RDA model (Borcard et al., 2011).

| RE SULTS
The mean foraging distance per 4 km 2 grid cell was 976 ± 268 (SE) meters, while the mean effective number of colonies per grid cell 12.7% of the variation was shared among the three categories.
In addition to our RDA models, we also carried out a series of multiple linear regressions to explore the combination of environmental parameters associated with the two genetic parameters representing high-quality bumblebee habitat. The effective number of colonies per grid cell was best predicted by a model (r 2 = 17.6%, Table 3) with the following 6 explanatory variables: relative area of aesthetic green space (grass/shrubs) (variance explained = 5.05%), relative area of roads (variance explained = 4.54%), "other" paved surfaces (variance explained = 4.36%), and relative area of watershed forests (variance explained = 3.38%) were all negatively associated with colNe, while city parks (variance explained = 7.09%) and relative area of watershed meadows (variance explained = 3.46%) all exhibited a positive relationship with colNe (Table 3). The log of mean foraging distance within a grid cell was best predicted by a model (r 2 = 35.0%) with four explanatory variables (Table 4): relative area of roads (variance explained = 15.5%), bare earth percentage (variance explained = 5.95%), and "other" paved surfaces (variance explained = 0.11%) were positively associated with foraging distance, while house density (variance explained = 18.8%) was negatively associated with foraging distance (Table 4).

| DISCUSS ION
Our high-resolution analyses of the landscape features that enhance pollinator habitats in urban environments in Toronto are consistent with previous studies on bumblebees and other generalist native bees. Specifically, previous studies linked impervious paved surfaces and availability of forage as key negative and positive drivers, respectively, of habitat quality for native bee communities (Bennett & Lovell, 2019;Birdshire et al., 2020;Carvell et al., 2011;Egerer et al., 2020;Fortel et al., 2014;Geslin et al., 2016;Goulson et al., 2010;Jha & Kremen, 2013a,b).
In our correlation analysis, the effective number of colonies (colNe) within a grid cell was negatively associated with the relative area of "other paved surfaces" (Figure 3). Both the RDA analysis ( Figure 4) and the multiple linear regression analyses (Table 3) also showed the benefits of green space and the disadvantages of man-made structures on colony density. These results are intuitive given that B. impatiens primarily nest in underground cavities that are often near trees and woody shrubs (Lanterman et al., 2019).

F I G U R E 3
Correlation structure of environmental, demographic, and population parameters. Only statistically significant correlations are depicted (after multiple testing correction), and the size and color of the dots represent the value of the Spearman correlation coefficient. Like colony density, the average foraging distance was also impacted by the degree of urbanization within Toronto. Our correlation analysis showed that more buildings and roads were associated with greater foraging distances (Figure 3). The RDA (Figure 4) and multiple linear regression analyses (Table 4) supported and expanded on the simpler analysis. In our RDA model (Figure 4), bumblebee foraging distance increased with higher densities of buildings, roads, paved surfaces, bare earth, and human population size (i.e., greater urbanization). Grid cells with more forests, city parks, and residential housing were associated with lower foraging distance. Our best fitting multiple linear regression model (Table 4) recapitulated the negative influence of urbanization on foraging distance: Both the relative area of roads, other paved surfaces, and bare earth were associated with higher foraging distances. Our findings are consistent with previous research showing that human structures, such as roads and railroads, substantially restrict foraging by bumblebees workers who avoid flying across these man-made structures (Bhattacharya et al., 2003). Surprisingly, the relative density of houses within grid cells was associated with shorter foraging distances. This likely reflects that the type of urbanization matters. In other words, in the absence of city parks or watershed forests, urban areas with a higher density of single or multiple family houses likely provide better foraging opportunities for bumblebees relative to urban areas with a high density of multilevel buildings. Houses in Toronto typically have a front yard and a back yard, which can provide some foraging opportunities for bumblebees throughout their active flight season (e.g., flowering trees, small-scale gardens, and weeds).
There has been conflicting evidence on the role of the "luxury effect" (Leong et al., 2018) in pollinator conservation. One study in the city of Chicago (USA) showed that lower income neighborhoods faired similarly in terms of bee abundance and diversity relative to higher income neighborhoods (Lowenstein et al., 2014). However, a F I G U R E 4 A redundancy analysis (RDA) for bumblebee foraging distance and effective number of colonies explained 55.2% of variation with an adjusted r 2 of 30.8% based on explanatory variables accounting for natural habitat, human-made infrastructure, and human population and demographic factors. RDA axis 1 explains 33.6% of the variation and RDA axis 2 explains 16.7% of the variation. The arrows in blue represent the explanatory variables, and the variables in red are the response variables. Notes: Grass/Shrub % has the same direction and weight as Elevation. Average foraging distance was log-transformed before analysis (see methods). Table 2 contains an appendix of terms   (Jha & Kremen, 2013a,b) and native bee communities (Bennett & Lovell, 2019;Birdshire et al., 2020;Egerer et al., 2020;Fortel et al., 2014;Geslin et al., 2016). Similarly, functional green space in Toronto was associated with shorter foraging distances in our RDA analysis, which is consistent with previous research on the importance of forage for many other native bumblebees (Carvell et al., 2011(Carvell et al., , 2012(Carvell et al., , 2017Goulson et al., 2010;Redhead et al., 2016). We thus think it is reasonable to assume that our analysis on B. impatiens can help us identify high-quality habitat for other native pollinators given that: (1) landscape features associated with high-quality habitat for B. impatiens in Toronto have also been associated with high-quality habitat for native bee communities elsewhere; (2) B. impatiens workers visit a diverse assemblage of plants that attract and provide resources to wide assemblages of native bees.
Overall, our study shows the importance of functional green space in providing high-quality habitat for bumblebees in an urban environment such as the city of Toronto. Both colony density and foraging distance were influenced by the degree of urbanization in Toronto. Impervious paved surfaces (e.g., roads, buildings, and other paved surfaces) were associated with differences in colony density and foraging distances. While green space is important for bumblebees, our study indicated that natural/functional green spaces, such as city parks and forests, were often beneficial relative to aesthetic green spaces, such as lawns. Our analysis suggests two simple strategies for improving bumblebee habitat within cities. First, conversion of paved surfaces to functional green space such as parks and meadows is likely to have a significant influence on the quality of pollinator habitats in Toronto. Second, our RDA ( Figure 4) suggests that converting aesthetic green space (i.e., lawns, which is orthogonal to the foraging distance vector and is opposite to the colony density vector) into more functional natural green space (e.g., flowering meadows, which tend to be associated with lower foraging distances) can improve the foraging opportunities of bumblebee colonies in Toronto.

ACK N OWLED G M ENTS
We thank Kailey Michnal for collecting bumblebee samples. This study was funded through a NSERC Discovery Grant to AZ, with support from York University's Research Chair program to AZ, SRC, and SS. We thank York University's Center for Bee Ecology, Evolution and Conservation for facilitating our collaborative research.

CO N FLI C T O F I NTE R E S T
None declared.

DATA AVA I L A B I L I T Y S TAT E M E N T
All raw and processed data can be found as Data S1. The code for carrying out the statistical landscape analyses can be downloaded on Github (https://github.com/mimri t/Toron to_Bumble_Bees).